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Abstract: We present an algorithm for the integrand-level reduction of multi-loop am- 
plitudes of renormalizable field theories, based on computational algebraic geometry. This 
algorithm uses (1) the Grobner basis method to determine the basis for integrand-level 
reduction, (2) the primary decomposition of an ideal to classify all inequivalent solutions 
of unitarity cuts. The resulting basis and cut solutions can be used to reconstruct the 
integrand from unitarity cuts, via polynomial fitting techniques. The basis determina- 
tion part of the algorithm has been implemented in the Mathematica package, BasisDet. 
The primary decomposition part can be readily carried out by algebraic geometry soft- 
wares, with the output of the package BasisDet. The algorithm works in both D = 4 and 
D = 4 — 2e dimensions, and we present some two and three-loop examples of applications 
of this algorithm. 
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1 Introduction 

The study of higher-loop amplitudes for gauge theories is important for both theoreti- 
cal and phenomenological reasons. The data analysis of Large Hadron Collider (LHC) 
physics requires great accuracy of the standard-model cross sections computation. For 
many channels, not only the next-to-leading order (NLO) amplitudes, but also the next- 
to-next-to-leading order (NNLO) amplitudes are important in order to control theoretical 
uncertainties. 

The traditional Feynman diagram approach for amplitude calculation becomes very 
complicated in the higher- loop cases. Integration- by-parts (IBP) identities were used to 
reduce the number of integrals in loop diagrams [1], and efficiently implemented in Laporta 
algorithm [2]. The method of Grobner basis was used to express Feynman integrals as a 
linear combination of master integrals [3-5] . 
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New methods based on unitarity [6-8] decompose loop amplitudes as the product of 
tree amplitudes and greatly simplify the computation. There is a particularly convenient 
method: integrand-level reduction (OPP method) [9], which decomposes the amplitude 
directly at the integrand level. OPP method can be used to automatically and efficiently 
calculate one-loop amplitudes with multiple legs [10-15]. The original OPP method can 
be generalized to D = 4 — 2e at one loop [16-18]. The polynomial division method was 
first used in the integrand-level reduction for one-loop amplitudes in ref. [19]. 

The generalized unitarity method also applies to two- loop amplitudes: Buchbinder 
and Cachazo calculated two-loop planar amplitudes for N = 4 super- Yang- Mills [20], by 
generalized unitarity. Gluza, Kajda and Kosower [21] used a Grobner basis to find IBP 
relations without double propagators, and then determined the master integrals for two- 
loop planar diagrams. The IBP relations can also be generated by linear algebra, without 
using Grobner basis [22]. With these master integrals, Kosower and Larsen [23] applied 
maximal unitarity method for two-loop planar diagrams to obtain the coefficients of master 
integrals from the products of tree amplitudes. Furthermore, Larsen [24] applied this 
method to two-loop six-point amplitudes, using multidimensional contour integrals. The 
two-loop double-box diagram maximal-cut solutions can be related to Riemann surfaces, 
whose geometry uniquely defines the contour integrals [25]. 

Alternatively, using integrand-level reduction, Mastrolia and Ossola [26] applied OPP- 
like methods to study two- loop N = 4 super- Yang-Mills amplitudes. Badger, Frellesvig and 
Zhang [27] then used the Gram-matrix method to find the integrand basis systematically 
for two- loop amplitudes of general renormalizable theories. They calculated the double- 
box and crossed-box contributions to two-loop four-point N = 0, 1, 2, 4 (super)-Yang-Mills 
amplitudes. 

It is interesting to generalize and automate the integrand-level reduction to amplitudes 
with more legs and more than two loops. The main limitations in the previous integrand- 
level reductions are, 

1. The basis for integrand-level reduction grows quickly as the number of loops in- 
creases. At three-loop order and the beyond, the Gram-matrix method becomes very 
complicated and the integrand-level basis is hard to obtain. 

2. It is difficult to find and classify all inequivalent unitarity cut solutions for compli- 
cated diagrams. It is necessary to find all cut solutions, to reconstruct the integrand. 
However, for diagrams with many legs or more than two loops, the solutions be- 
come complicated. And often two solutions appear to be different, are but actually 
equivalent after reparametrization [25]. We have to remove this redundancy before 
reconstructing the integrand. 

These difficulties come from the complexity of the algebraic system of cut equations. The 
ideal approach to deal with these problems is computational algebraic geometry. In this 
paper, we reformulate these two problems as classic mathematical problems and solve them 
by powerful mathematical tools, 
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1. Integrand- level basis is equivalent to the linear basis in the quotient ring, of poly- 
nomials in irreducible scalar products (ISPs) modulo the cut equations. Then the 
integrand basis can be generated automatically using the standard Grobner basis and 
polynomial reduction methods [28, Ch. 5]. 

2. The collection of all cut solutions is an algebraic set. The latter can be uniquely de- 
composed to a finite number of affine varieties. Each variety is an independent solu- 
tion of the unitarity cuts, and different varieties are not equivalent by reparametriza- 
tion. In practice, this decomposition is automatically done by primary decomposition 
of an ideal [29, Ch. 1]. This classifies all inequivalent unitarity cut solutions. Further- 
more, dimension theory in algebraic geometry [29, Ch. 1] can determine the number 
of free parameters in each solution. 

We implement the first part of our algorithm in the Mathematica package BasisDet 
which can automatically generate the integrand-level basis. It also provides a list of ir- 
reducible scalar products (ISP)'s and the ideal / generated by the cut equations. The 
latter information can be directly used by computational algebraic geometry software, like 
Macaulay2 [30] , to carry out the second part of the algorithm. Once the primary decom- 
position is done, we get all inequivalent solutions of the unitarity cuts. Furthermore, for 
each solution, the software will find the number of free parameters. 

The package BasisDet has been tested for D = 4 and D = 4 — 2e one-loop box, trian- 
gle and bubble diagrams, D = 4 two-loop four-point double-box, crossed-box, pentagon- 
triangle diagrams, D = 4 two-loop five-point double-box diagram, pentagon-box diagram, 
and D = 4 — 2e two-loop four-point diagram. It has also been tested in two-loop level 
diagrams beyond maximal unitarity, for example, D = 4 two-loop four-point box-triangle, 
sunset and double-bubble diagrams. The output bases have been verified for all these cases. 

We have also used this algorithm to calculate D = 4 three-loop triple-box basis, and 
have verified that terms inside the basis are linearly-independent on the unitarity cuts. 
We also successfully carried out a primary decomposition on this diagram to find all the 
inequivalent cut solutions. 

This paper is organized as follows. In section 2, we briefly review the known integrand- 
level reduction for one and two-loop diagrams. The limitation of previous approaches is also 
pointed out. In section 3, our new algorithm is presented, and its validity is mathematically 
proven. Then in section 4, several examples are presented for one, two and three-loop 
diagrams. Finally, our conclusion and discussion on future directions are provided in section 
5. The manual for the package BasisDet is given in Appendix A. 

The package BasisDet and examples are included in ancillary files of the arXiv version 
of this paper. The package and its future updates can also be downloaded from the website, 
http : / /www . nbi . dk/ ~zhang/ BasisDet . html. 

2 Review of integrand-level reduction methods 

In this section, we briefly review integrand-level reduction for one and two-loop amplitudes. 
(See [18] for detailed review of the one-loop integrand reduction.) 
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2.1 One-loop integrand- level reduction 
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Figure 1. One-loop box diagram 
Schematically, for D = 4, an one-loop amplitude must be decomposed as [9], 
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(2.1) 



where we define the propagators Z)^ = (fe— Pii,i x -i) , Pi,j = Ylk=iPk (figure 1) as the sum of 
external momenta such that p^^-i = 0, and have taken the restriction that all propagators 
are massless. We must require that A4 i j 1 j 2 j 3 j 4 contain no term which is proportional to Di x , 
x = 1, ... 4, otherwise one of the denominator in the integral is cancelled out. Similarly, 
A3 j i 1 j 2 i 3 must contain no term proportional to D^, x = 1, ... 3 and so on. 

Consider A4 i j 1 j 2 j 3 j 4 first. There exists a vector u) perpendicular to p^^—i, x = 2,3,4. 
Of all the scalar products involving loop momenta, only k ■ oj is not a polynomial in de- 
nominators Djj, Di 2 , Di 3 and Di A . We call such scalar products irreducible scalar products 
(ISPs) and the other scalars products reducible scalar products (RSPs). A^i 1 i 2 i 3 i i should 
be a function of ISPs, i.e. (k ■ oj) only. 



Furthermore, we find that A4 j j 1 j 2 j 3 j 4 is at most linear in k ■ 

&i,ilhhi A = CO + Ci(fc • Oj). 



(2.2) 



Co and c\ are constants independent of the loop momentum. Higher-order terms in (k ■ oj) 
are absent, because 



(k ■ oj) 2 = linear combination of Di x , D{ 2 , Di z and Di 



(2.3) 



- 4 - 



The coefficients cq and c\ can be calculated from quadruple cuts, Di 1 (k^) = Di 2 (k^) = 
D i3 (kW) = D i4 (k^) = 0, where s = 1,2 since there are two cut solutions. The integrand 
at the two cut solutions determined the two coefficients Co and c\. Note that although the 
term c\{k ■ uj) integrates to zero, it is necessary to keep it for the triple-cut calculation. We 
call the set {l,(k ■ u)} the integrand basis for D = 4 one- loop quadruple cut and {k ■ uj) 
the spurious term. 

Similarly, A.^^^ can be reconstructed from triple cuts. We have two vectors uj\ and 
0J2, which are in perpendicular to the external momenta. There are two ISPs, k ■ uj\ and 
k ■ uj 2 . The expansion over integrand basis reads, 

A 3,ui2i 3 ( fc ) = c oo + cio(& • wi) + coi(fc • uj 2 ) + c u (k ■ ui)(k ■ u) 2 ) 
+c 12 (k ■ wi)(/c • lo 2 ) 2 + c 2 \{k ■ ui) 2 (k ■ u 2 ) + c 20; 02 ((k ■ uji) 2 - (k ■ u 2 ) 2 ) . (2.4) 

The basis contains 7 terms, of which 6 are spurious. There are two cut solutions for the 
triple cut, 

Ax(fc W (r)) = D l2 {k^( T )) = D i3 (k^(r)) = 0, (2.5) 

with s = 1,2, but each of them contains one complex free parameter r. The original 
numerator of the integral at triple cuts, with all A^j^^ subtracted, becomes a Laurent 
series in r. The corresponding Laurent coefficients determine all the 7 "c" coefficients of 
eq. (2.4). 

The one-loop integrand-level reduction also works for D = 4 — 2e [18]. The loop 
momenta contains both the four-dimensional part and the extra-dimensional part , 

k = k^+k ± , {k L ) 2 = -ti 2 . (2.6) 

For the quadruple cut, the basis has larger size than that of the D = 4 case. Instead of 
(2.3), we have 

(k ■ uj) 2 — /x 2 = linear combination of D^, Di 2 , Di 3 and D{ 4 . (2.7) 

So we can remove either (k ■ uj) 2 or fj? to obtain an integrand basis. One convenient choice 
is 

A tiili 3 u = + ci{k ■ uj) + c 2/ u 2 + c 3 {k ■ uj)fi 2 + c 4 /A (2.8) 

The D = 4 — 2e quadruple cut has only one solution. This solution depends on one complex 
free parameter r. Note that geometrically, this solution is complex one-dimensional, and 
contains the two D = 4 box quadruple cut solutions (zero-dimensional) as two isolated 
points. The Taylor expansion in r of the integrand, at the quadruple cut, determined the 
coefficients Co, c±, c 2 , C3 and C4. 

throughout this paper, we use the scheme that all external momenta and polarization vectors have no 
(— 2e)-dimensional components. 
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2.2 Two-loop integrand-level reduction 

For the one-loop cases considered above, it is relatively easy to determine the integrand 
basis and find the unitarity cut solutions. However, in the two-loop cases, it is much harder 
to find the integrand basis and the unitarity cut solutions are more complicated. 

Mastrolia and Ossola [26] applied the OPP-like method for two- loop N = 4 super Yang- 
Mills amplitudes. Two-loop four-point amplitudes for general renormalizable theories were 
calculated in [27]. To show clearly new features of two-loop integrand-level reduction, we 
review the Gram-matrix method presented in ref. [27]. 

For example, consider the two-loop four-point planar diagram (Figure. 2). The integrand- 
level reduction reads, 
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(2.9) 



where . . . stands for the integrals with less than 7 propagators. Our aim is to recon- 
struct the double box function A^u*3A* from hepta-cuts (maximal cut for diagrams with 
7 propagators). Again, there exists one vector uj which is perpendicular to all the external 
momenta. There are four ISPs: (k ■ p^), (q ■ p±), (k ■ cj) and (q ■ uj). The integrand basis 




Figure 2. Four-point two- loop planar diagram 

consists of terms with the form, 

{k-ptTiq-pyTik-uiYiq-uY, (2.10) 

where to, n, a and /3 are non-negative integers. For renormalizable theories, power counting 
requires that, m + n + a + (3<6, m + a < 4 and n + /3 < 4. Furthermore, it is easy to see 
that (k ■ to) 2 , (q ■ to) 2 and (k ■ u))(q ■ ui) are linear combinations of the seven denominators. 
Hence a < 1, /3 < 1 and a ■ (3 = 0. The above analysis is similar to that of one- loop cases, 
and the integrand basis appears to contain 56 terms. 

However, there are more constraints. For four dimension momenta, the determinants 
of 5 x 5 Gram matrices are zero, 

, „ ( 12 4 k q\ n , /l24£:<A n , ^ f 1 2 4 k q\ n /n 
detG =0, detG =0, detG y = 0. (2.11) 

yi24kqj yi24ukj \1 2 4 u q J 
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For example, on the hepta-cut, the first Gram-matrix relation reads, 

= 4(fc • p 4 ) 2 (q ■ Pl f + 2s(k ■ p 4 ) 2 (q ■ Pl ) + 2s(k ■ p 4 )(q ■ Pl ) 2 - st{k ■ p 4 ) (q ■ Pl ) (2.12) 

which means that either m < 2 or n < 2. Finally, the integrand basis for the double box 
amplitude is, 

A^° 2 X *34*(M) = E Cmn(a + 2p)(k-p4) m (q-piT(k-LUr(q-uf. (2.13) 

mna/3 

There are 16 non-spurious terms, 

(C000! CoiO, Cioo, C020i Cuo, C200) C030> C±20, C210 > C30O ' C 040 > C130 , C310, C400, Cuo, C410), (2-14) 

and 16 spurious terms, 

(cooi> con, cioi, cm, C201, C211, C301, C311, C002, C012, C102, C022, C112, C032, C122, C132). (2.15) 

We comment that, alternatively, (2.12) can be obtained using the elimination method on 
the 7 cut equations. However, its computation is quite long and not systematic, comparing 
with the Gram-matrix method. In two-loop cases, the Gram-matrix method provides a 
very efficient way to determine the basis. 

Again, we can determine the 32 coefficients from the hepta-cuts solutions. The solu- 
tions are much more complicated than one- loop cut solutions: there are 6 solutions, and 
each of which depends on a free parameter, r. From the Taylor or Laurent expansion of 
the integrand at hepta-cuts, we can solve for the 32 coefficients, in cases of Yang-Mills 
and TV = 1,2,4 super- Yang-Mills theories. Then IBP relations can further reduce the 16 
non-spurious integrals to two master integrals. However, to get the lower cut functions 
like the hexa-cut case, we have to subtract all the 32 terms first, not only the two master 
integrals. 

The four-point non-planar function A^j^^ has been determined by the same method 

[27]. 

The Gram-matrix method becomes more complicated as we attempt to add more loops 
and legs. Furthermore, it is not easy to automate the Gram-matrix method: once a Gram- 
matrix relation is found, we need to determine which monomial inside the relation should 
be removed from the basis. For diagrams with many legs or more than two loops, it is also 
complicated to classify all the cut solutions. Hence a new automatic algorithm is needed, 
to carry out integrand- level reduction for higher-loop and many- leg amplitudes. 

3 The algorithm 

We present an automatic algorithm for integrand-level basis determination for generalized 
unitarity, based on the techniques of computational algebraic geometry. The goal is (1) to 
find the integrand basis by Grobner basis method (2) to classify all inequivalent unitary 
cut solutions by primary decomposition and find the dimension of each solution. 
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3.1 Setup 

We parametrize the loop momenta using scalar products. This is a variation of van Neerven- 
Vermaseren basis [31]. This parameterization has the following advantages: 

• It does not depend on spinor helicity formalism or particular basis choices. 

• The cut equations in terms of scalar products have a particularly simple form. This 
makes it convenient to carry out primary decomposition later. It is also easier to 
apply polynomial fitting techniques to reconstruct the integrand. 

Consider an L-loop n-point diagram. The dimension is D = d or D = d — 2e. d 
is an integer which stands for the dimension of the physical spacetime, while — 2e is the 
number of extra dimensions introduced by dimensional regularization. In most examples, 
we consider d = 4. 

Let be the loop momenta and (pi,...p n ) be the external momenta. We 

use the scheme that all extra momenta and polarization vectors have no extra-dimensional 
components. The momenta pj can be either massless or massive. 

We choose a basis {ei, e^} for the physical spacetime. Each is either an external 
momentum or an u)j, that is a momentum perpendicular to all the external legs. We define 
the d x d Gram matrix, 

G d = G\ e i'-'**y (3.1) 

Gd is nonsingular, as it should be. 

For the case D = d — 2e, we decompose the loop momenta into physical and extra- 
dimensional components, 

k = l? + lt (3.2) 

and we define /iy = — if- ■ For the case D = d, we simply set if- = and all the ynj are 
absent. 

We parametrize l f^ using scalar products, (k ■ ej), 1 < j < D, 



(3.3) 



l? = (e 1 ,...,e d )Gj 1 \ 

\ (k ■ e d ) 

We define the set of (fundamental-)scalar products (SPs) to be 

SP = {(li • &j)|l < i < L,l < j < d}U {/%|1 < i < L,i < j < L}, D = d - 2e, (3.4) 

or 

SP = {(k ■ ej)\l < i < L,l < j < d}, D = d. (3.5) 

All the other scalar products, like (li ■ u), if, (li ■ lj), where u is a constant vector in 
the physical dimension, can be written as polynomial functions of (fundamental-)scalar 



-8- 



products, using the Gram matrix G^- For example, 

/(u- ei )\ 
k ■ u = {(k ■ ei), (k ■ e d )) G^ 1 : 

( (h • e i 



(3.6) 



k -lj = {{k ■e 1 ),...,(l i •e d ))G f2 1 



Hij. (3.7) 



Next, we consider the m-fold unitarity cut of the amplitude, i.e., m propagators are 
set to zero. 

D 1 (h,...l L ) = ... = D m (h, ...l L ) = 0, (3.8) 

Using the Gram matrix Gd, these cut equations can be expressed m polynomial equations 
in the SPs. We denote the polynomial ring of SPs, i.e. the collection of all polynoimals in 
SPs, as R'. Then we introduce the concept of an ideal in a ring [28, Ch. 1]: in general, an 
ideal J generated by several polynomials /i, ... in a ring S, is the subset of S, 

J =(/!,.. . f k ) = {oi/i + . . . + a k f k \V ai g S, 1 < i < k}, (3.9) 

where a^s are arbitrary polynomials in S. Here we define, 

/' = (£>!,... D m ), (3.10) 

which is the ideal generated by all the cut equations in terms of SPs. 

Some scalar products' values are uniquely determined at all cut solutions, i.e., they are 
polynomials of propagators, 

x = const + (£>!,..., D m ). (3.11) 

We may call these scalar products reducible scalar products (RSPs) and and all the other 
scalar products in §P irreducible scalar products (ISPs). 

In practice, we can extend the definition of RSPs. For example, if two scalar products 
x\ and 22 1 satisfy the relation, 

a\x\ + a 2 x 2 = const + 0(D 1 , . . . , D m ), (3.12) 

where a\ and a 2 are nonzero constant. We may pick up one of the two scalar products as 
RSP, say x±, and write it as a linear function of x 2 on the multiple cut. 

Hence we have the following formal definition of reducible scalar products (RSP) and 
and all the other scalar products irreducible scalar products (ISP): 

Definition 1. The set ISP of irreducible scalar products is a minimal subset o/SP, such that 
all the other scalar products can be expressed as linear functions in ISPs on the unitarity 
cut. 
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This definition minimizes the number of ISPs, so the following calculation will be 
simpler. The choice of ISP is not unique but different choices are equivalent. We have the 
following decomposition: 

SP = RSPUISP. (3.13) 

To simplify notations, we label the ISPs by x±, . . . x ni . 

We can eliminate all the RSPs from the cut equations to obtain a new set of algebraic 
equations F in ISPs. With an abuse of notations, 

F = {D k ( Xl , . . . x ni ) = 0|1 < k < m}, (3.14) 

where is the polynomial in ISPs obtained from rewriting the /c-th propagators in terms 
of ISPs, after RSPs are eliminated. 

We denote the polynomial ring of ISPs as R, and the ideal generated by D&S as /, 

I = (D 1 ,...,D m ) C R, (3.15) 

where (...) stands for an ideal generated by several polynomials. 

It is easy to identify the RSPs and ISPs by hand for one and two-loop diagrams. 
However, this calculation becomes messy for more complicated diagrams. In practice, the 
identification of the RSP and ISP can be done quickly and systematically using Grobner 
basis method, as described in appendix B. 

The algebraic equation system F in ISPs plays the central role in our algorithm. We 
will see that it contains all the information on cut solutions and the integrand basis. 

3.2 Algorithm for integrand basis determination 

In this section, we present an automatable algorithm for the determination of the integrand 
basis. 

From the previous section, we see that all Lorentz invariants can be reduced to poly- 
nomials of (fundamental-)scalar products. Furthermore, RSPs can be reduced to constants 
or linear functions of ISPs. Hence, schematically, on m-fold unitarity cuts of a L-loop 
amplitude, the numerator of the integrand is reduced to a polynomial of ISPs, like (2.2), 
(2.4), (2.8) and (2.13), 

^m l0 ° P = 5^ c Ql ... an ^3; 1 1 . . . x ni ' . (3.16) 

(ai,...,a„j)£S 

Here the tuple (ai, . . . a ni ) groups together the non-negative integer powers of the ISPs. 
We further require that Am loop have no dependence in the propagators Dj, i = 1, . . . , m. In 
previous examples, this was achieved by using cut equations directly for one-loop topolo- 
gies, or Gram-matrix method for two-loop topologies. The finite set S contains all the 
power tuples for the reduction. The constant coefficients c ai ... an are independent of loop 
momenta. They can be fitted from tree amplitudes by unitarity as in previous examples. 

The set B of the monomials x" 1 . . . x"™ 7 appearing in (3.16), is called the integrand 
basis. The goal is to determine this basis, or equivalently, the finite set S. We translate 
the requirement that Am' oop have no dependence in the propagators Di, i = 1, . . . , m into 
mathematical language, 
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Proposition 1. The monomials in the integrand basis must be linearly independent in the 
quotient ring R/I. 

Proof. Otherwise, there exist constant coefficients d ai ,,, ani , (a\, . . . ,a ni ) £ S such that 

m 

^2 dai.-.a^Xi 1 ■ ..Xn7 = ^fkD k (3.17) 

(ai,...,Q ni )eS k=l 

where fk's are polynomials in ISPs. Suppose that one coefficient dp 1 ...p ni is not zero. We 
define a subset S = S — {(/3i . . . f3 ni )}. Then Am' oop can be reduced even further, 



a L- loop _ \ ^ctx 

(ai,...,a nj )€S 



c ai ... a - d ai ...a ni )x l 1 ...x n / + > f k D k (3.18) 

Thus Am" loop still depends on D\, . . . D m . We can redefine the first term in (3.18) as 
Am loop , § as the new power set for the basis, and move the second term in (3.18) to 
fewer-propagator integrals. The size of the basis decreases by one, so the reduction is not 
complete. □ 

There is a classic method to find the complete linearly independent basis in R/I: 
Buchberger's algorithm [28, Ch. 5]. (See [28, Ch. 2] for a review of Grobner basis.) 

1. Define a monomial ordering in R and calculate the corresponding Grobner basis G(I) 
of /. Denote LT(K) as the collection of leading terms of all polynomials in a set K, 
according to this monomial ordering. 

2. Compute LT(G(I)), the leading terms of all polynomials in G(I). Obtain (LT(G(I))) , 
the ideal generated by LT(G(I)). By the properties of Grobner basis, (LT(I)) = 
(LT(G(I))) , where (LT(I)) is the ideal generated by all leading terms in I. 

3. Then the linear basis of R/I is B, which is the set of all monomials which are not in 
(LT(I)). 

B = {x?...x% I \x?...x%> ?{LT(I))}. (3.19) 

This method is fast. However, the basis generated by Buchberger's method usually 
contains an infinite number of terms, since the renormalizablity conditions have not been 
imposed. We find that after the ring R is reduced to B, it is not easy to impose renormal- 
izablity conditions. 

Hence, we propose the following alternative algorithm for basis determination, based 
on multivariate synthetic division, 

1. Define a monomial ordering in R and calculate the corresponding Grobner basis G(I) 
oil. 
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2. Generate the set A of all monomials in ISPs, which satisfy renormalizablity condi- 
tions. A must be a finite set. 

3. For each monomial a,j(x\ . . . x ni ) in A, 1 < j < \A\, 

• Carry out the multivariate synthetic division of a,- by the Grobner basis G{I). 

aj( Xl .. 

• x ni) — 9j 

(xi...x ni ) +Vj(xi.. 

■ %nj ) j 9j 

( Xl ...x ni )el (3.20) 

where rj(x± . . .x n[ ) is the remainder of multivariate synthetic division. Given 
the Grobner basis G(I), rj(xi . . . x ni ) is uniquely determined. 

• Decompose Vj(x\ . . . x ni ) as monomials and collect them in a set Bj. 

4. The integrand basis B is then, 



The validity of this algorithm can be verified as follows, 

• The monomials in B are linearly independent in R/I. Multivariate synthetic division 
by Grobner basis ensures that all monomials in Bj are not in (LT(I)), therefore B 
is a subset of B. So by a corollary of Buchberger's method, linear independence is 
proven. 

• The basis B is big enough for integrand-level reduction. From step 3, we see that 
every renormalizable term in the numerator of the integrand is reduced to monomials 
in B. In other words, it is a sum of a linear combination of monomials in B and other 
terms vanishing on the unitarity cut. 

We implement this part of our algorithm in the Mathematica package BasisDet. The 
Grobner basis calculation and multivariate synthetic division are done by the functions in 
Mathematica. The monomial order is chosen as "degree lexicographic" ( "deglex" in math- 
ematica language, see [28, Ch. 2] for a review of monomial orderings.) and the coefficient 
field can be chosen as rational functions for analytic computation, or rational numbers for 
numeric computation. 

3.3 Primary decomposition of cut solutions 

Given the cut equations in ISP variables, or equivalently, the ideal /, the following questions 
naturally arise: 

• How many inequivalent cut solutions are there? 

• For each cut solution, how many free parameters are needed to parametrize it? In 
the other world, which is the dimension of each cut-solution? 




(3.21) 



j 
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These questions can be studied systematically using algebraic geometry. We again translate 
these problems to mathematical language. Consider the affine space A ni = (xi, . . . x ni ). 
The ideal I defines an affine algebraic set Y in A Hl , [29, Ch. 1] 

Y = Z(I) = {(z 1 ,...z ni )\D k (z 1 ,...z ni ) = 0, V/c} (3.22) 

which is the collection of all cut solutions in term of ISPs. 

In general, Y can always be decomposed uniquely to the union of a finite number of 
irreducible components [29, Ch. 1], 

"sol 

Y = (J Y a , Y a (/L Y b , if a + b. (3.23) 

a=l 

where each Y a is an affine variety. Here, n so \ is the number of irreducible components. Dif- 
ferent components are not related by parameter redefinition. Each irreducible component 
corresponds to a cut solution. So we have the following proposition: 

Proposition 2. The inequivalent cut solutions are in one-to-one correspondence with the 
irreducible components of the algebraic set Y. In particular, the number of inequivalent cut 
solutions equals the number of the irreducible components of the algebraic set Y. 

This decomposition can be achieved easily by the algebraic method, primary decom- 
position of an ideal [29, Ch. 1]. Since R (the polynomial ring of ISPs) is a Noetherian ring, 
the primary decomposition of I uniquely exists (Lasker-Noether theorem) [32], 



Pi la- (3.24) 



o=l 



where s is a finite integer and each I a is a primary ideal. Furthermore, the primary 
decomposition guaranteed that, 



'Ia^Vh, if a^b. (3.25) 



where y[T a is the radical of I a . 2 Because I a is primary, ~JT a is a prime ideal. 

Hence we have the corresponding decomposition of Y . Define Z(I a ) to be the zero-locus 
(set of all solution points) of the ideal I a , 



Y = Z(I) = (J Z(I a ) = |J Z{^T a ). (3.27) 



a=l a=l 



Since \fT~ a is prime, Z(I a ) = Z(\/ r T2) is an affine variety [29, Ch. 1] which is irreducible. 
Then we define n so i = s and Y a = Z{yfl^), and the decomposition is done. 



2 The radical of an ideal J is the set of all elements a, such that a n £ J, where n is some positive integer. 
[28, Ch. 4] 
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The dimension of each component is given by dimension theory in commutative algebra 
[29, Ch. 1], 

dim Y a =nj- height (V^o)- (3-28) 

where height ("/la) is the height of the prime ideal y/Ta~, which is defined to be the largest 
integer N, for all possible series of prime ideals (0) = po C p\ C . . . C pn = \fla- Recall 
that SP is the set of (fundamental-)scalar products, which is defined in (3.4). Note the 
dimY^ may not equal |SP| —m, the difference between the number of (fundamental-)scalar 
products and the number of cut equations, because of the possible redundancy in the cut 
equations. Furthermore, for a ^ b, dim Y a may not equal dim Yj,, since they are independent 
components. 

Once all the irreducible components are obtained, we can parametrize each inequivalent 
solution. Together with the RSPs, the explicit form for loop momenta li at each cut solution 
can be recovered. 

The primary decomposition (3.24) and dimension (3.28) can be calculated using com- 
putational algebraic geometry software, for example, by standard built-in functions in 
Macaulay2 [30] by Daniel Grayson and Michael Stillman. Alternatively, if we only need 
the number of irreducible components, then a numeric algebraic geometry approach could 
be applied, as described in [33]. 

4 Examples 

We implemented the basis determination part of our algorithm in the Mathematica package 
BasisDet. The only required inputs are the kinematic relations for external legs, a list of 
propagators and the renormalization conditions. The output is the integrand basis. It 
also provides /, as in (3.15), the ideal generated by the cut equations in terms of ISPs. 
Then we can carry out the primary decomposition and dimension theory calculation in the 
computational algebraic geometry program, Macaulay2, with the ideal I obtained from 
BasisDet. Here we list several examples of application of our algorithm. All computations 
were done on a laptop with an Intel core i7 CPU. 

4.1 D = 4 — 2e one-loop four-point box topology 

Take D = 4 — 2e, and consider the one-loop contribution with box topology to four-point- 
all-massless amplitude. The BasisDet package takes 0.05 seconds to generate the basis in 
the analytic mode (see the appendix A for the modes of the package) , 

A ti3Ui4 = CO + Cl(k ■ + C 2/ U 2 + C 3 (k • U)fl 2 + C4/A (4.1) 

which is exactly the same basis as (2.8), which was obtained in ref. [18]. The package 
automatically find the two ISPs (k • uj) and \j? . The cut equations, after all RSP are 
eliminated, become one equation, 

4 (fc • u) 2 - —fi 2 -t 2 = 0, (4.2) 
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where s, t, u are Mandelstam variables. It is clear that this equation defines an irreducible 
parabola in the parameter space (k-u, fx 2 ). So there is only one solution, with the dimension 
1. As a trivial test, we can also see that from primary decomposition. Let, 

I=U(k-uf-— ii\-t 2 ). (4.3) 
s 

Macaulay2 determines that / itself is primary, so no decomposition is needed. It also 
automatically finds that dim 1=1, which means there is one free parameter for the solution. 

4.2 Two-loop examples 

First, we consider D = 4 four-massless-particle amplitude with two- loop double-box topol- 
ogy (Figure. 2). The BasisDet package takes 0.95 seconds to generate the basis in the 
analytic mode, or 0.43 seconds to generate the same basis in the numeric mode. 

AgWfc,g)= E Cmn( a+ 2p)(k-P4) m (q-pi) n (k.ur(q-oo)P. (4.4) 

mnaf} 

with 16 non-spurious terms, 

(COOCH CoiO, Cioo, C020, Cuo, C200; C030, Cl20i C210, C300, C()40> C130, C310, C400, C140, C410), (4.5) 

which are exactly the same as (2.14). There are also 16 spurious terms, 

(C001) COll) Cioi, C201, C30I, C002; C012, C022, C032, C102, C112, C122, C132, C202, C302, C402)- (4.6) 

Note that although the number of terms is the same as (2.15), some terms are different 
from (2.15). It means we get a different but equivalent integrand basis. Even for the one 
loop D = 4 — 2e box quadruple cut, there are already several different choices of basis. 
We can check explicitly that the difference between (4.6) and (2.15) is proportional to the 
seven propagators, so it does not change the double-box contribution to the amplitude. 

There are four ISPs, {l\ ■ p A ), (l 2 • pi), (h • oS) and (l 2 • w). The cut equations in ISPs 
read, 

h = ~t 2 + 4t(/i • PA ) - 4(/i • p A ) 2 + A(h ■ ^) 2 = 0, (4.7) 
f 2 = -t 2 + 4i(/ 2 • pi) - 4(Z 2 • Pl ) 2 + 4(Z 2 • u) 2 = 0, (4.8) 
f 3 = s (-(h • p A ) 2 - 2(Zi • p A )(l 2 ■ Pl ) + (h • u) 2 + 2{h ■ u){l 2 ■ u) - (h ■ Pi? + (h ■ u) 2 ) 

-4t(h ■ p 4 )(l 2 ■ Pl ) = 0. (4.9) 

In this case, the ideal I = (/1, /s) is quite complicated. It is not easy to find the inequiv- 
alent solutions by hand or by elementary analytic geometry. We use primary decomposition 
to find inequivalent solutions automatically, for example, in Macaulay2, in just a couple of 
seconds, 

6 

I=f\li, (4.10) 
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where IjS are six primary ideals: 



h = (h • Pi, 2{l 2 • w) - 2(/ 2 • Pl ) + t, 2(h .u)-t) (4.11) 

/ 2 = (h ■ p 4 , 2{l 2 ■ u) + 2(l 2 ■ pi) - t, 2(h -u) + t) (4.12) 

h = (h ■ Pi,2(h -u) + t, 2(h ■ u) + 2(h • p 4 ) - i) (4.13) 

h = (h ■ PiXh • u) - t, 2(h ■ u) - 2(h ■ p 4 ) + t) (4.14) 
I s = (2(/ 2 • w) - 2{l 2 ■ p x ) + t, 2(h ■ oj) - 2(h ■ p 4 ) + t, 

4(/i • p 4 )(^2 • Pi) + 2{h ■ p 4 )s + 2(/ 2 • Pl )s - st) (4.15) 
7 6 = (2(Z 2 • u) + 2{l 2 ■ px) - t, 2(h ■ oj) + 2(h ■ p 4 ) - t, 

4(Zi • p 4 )(/ 2 • pi) + 2(ii • p 4 )s + 2(Z 2 • Pl )s - st). (4.16) 



So there are 6 inequivalent unitarity cut solutions, consistent with [23, 26]. Furthermore, 
Macaulay2 automatically finds that every solution of Ii has dimension 1. 

Note that all ij's are generated by simple polynomials, so it is straightforward to solve 
them for ISPs. Then using the Gram matrix relation, (3.3), we can rewrite the solutions in 
terms of the loop momenta lx and l 2 and find the one-to-one correspondence with the six 
solutions in ref. [23]. However, this step is not necessary since we can fit the coefficients 
c mn(o+2/3) directly from the solution for ISPs, as described in ref. [27]. 

Similarly, we can apply the same method on other two-loop diagrams using BasisDet 
and Macaulay2. Several examples are listed in Table 1. 





Diagram 


#ISP 


n-NS 


n s 


^basis 


^Solution 


Box-triangle 


\ 
/ 


s 
/ 




5 


18 


51 


69 


4 


Five-point double-box 


\ 


/ 


4 


32 





32 


6 










/ 




\ 


Sunset 




8 


12 


30 


42 


1 


Double-bubble 






6 


8 


48 


56 


1 



Table 1. Several examples of the integrand- level reduction of D = 4 two- loop diagrams: All 
external legs are massless. "#ISP" is the number of ISPs, tins an d «s are the numbers of non- 
spurious and spurious terms in the integrand basis, respectively, nbasis = n-NS + n s is the total 
number of terms. "^Solution" is the number of inequivalent solutions. The explicit expression of 
the integrand basis can be obtained by running the code "example. nb" with BasisDet. 
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4.3 D = 4 three-loop triple-box topology 

Consider D = 4, four-massless-particle diagram with three-loop triple-box topology (figure 
3). The package uses about 42 seconds in numeric mode or 4 minutes in analytic mode, to 




Figure 3. Four-point three- loop planar diagram 

generate the same integrand basis. It contains 199 non-spurious terms and 199 spurious 
terms. 

Furthermore we used Macaulay2 to find the inequivalent cut solutions by primary 
decomposition. It takes about 2 minutes to get, 

14 

8=1 

so there are 14 inequivalent solutions. And, 

dim/i = 2, z = l,..., 14. (4.18) 

Every solution thus depends on two free parameters. These solutions have been verified 
both analytically and numerically. Furthermore, by the explicit solutions, we can check 
that the 398 = 199 + 199 terms in the basis are linearly independent on the unitarity cuts. 
This validates the basis. 

With the integrand basis and all inequivalent solutions, we can reconstruct the triple- 
box contribution to three-loop amplitude for any renormalizable theory, via the polynomial 
fitting techniques. 

5 Conclusions and future directions 

In the paper, we have presented a new method for integrand-level reduction, based on 
computational algebraic geometry. It applies (1) a Grobner basis to find the basis for 
integrand-level reduction, (2) a primary decomposition of ideals to classify all inequivalent 
solutions of unitarity cuts. The first part is realized in our Mathematica package BasisDet, 
which automatically generates the basis from the propagator information. This package 
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also generates the ideal of the cut equations, which can be used as the input of primary 
decomposition. Then computational algebraic geometry software, like Macaulay2 [30] can 
classify all inequivalent cut solutions and determine the dimension of each solution. 

Since this method has no dependence on the spacetime dimension, the number of loops 
or the number of external legs, it works for general multi-loop diagrams of renormalizable 
theories. 

We applied this method to many one-loop and two-loop topologies. We have also used 
it to generate the correct basis and cut solutions for the three-loop triple-box topology. This 
method presented in this paper can be used to calculate many two-loop and higher-loop 
amplitudes, via polynomial fitting techniques. 

In future, it would be interesting to work on the following directions, 

• Symmetries in the diagram. It is interesting to find a monomial ordering to keep 
symmetries in the diagram manifest, for the basis determination part of our algo- 
rithm. Then this algorithm could be sped up considerably if all symmetries could be 
made manifest. 

• Automatic parametrization of each cut solutions. We would like to find an auto- 
matic way of parametrizing each cut solution, after the dimension of each solution 
is obtained. It would be helpful for the polynomial fit process, to reconstruct the 
integrand from unitarity cuts. 
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A Manual for the package BasisDet 1.01 

This current code is powered by Mathematica with its embedded Grobner basis function. 

The latest version of the package is available on the website, http : / / www . nbi . dk/ ~zhang/ BasisDet . ht 

A.l Set up 

The main program is "BasisDet-a-b.m" , where a. b is the version number of the package. 
It should be executed as, 

«' ' /path/BasisDet-a-b.m' ' 

where "/path/" is the path for the package. 
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A. 2 Input for loop diagrams 

The following variables need to be denned for the basis determination, 

• L. It is the number of loops in the diagram. 

• Dim. It is the dimension of the spacetime, which should be d or d — 2e. d is a positive 
integer and in most cases d = 4. 

• n. It is the number of external legs. 

• ExternalMomentaBasis. It is a list of external momenta in the basis for physical 
spacetime. Note that for n-point amplitude, because of momentum conservation, we 
can pick up at most n — 1 external momenta for the basis. In summary, 

— If n < d + 1, we need to put n — 1 external momenta in "ExternalMomentaBa- 
sis" . The program will automatically name the d — n + 1 spurious vectors as 

i0\, . . . ,LUd-n+l- 

— If n > d + 1, we need to put d external momenta in "ExternalMomentaBasis". 

• Kinematics. This is the list for replacement rules, from the kinematics. Note that 
only the scalar products of vectors in the basis need to be defined. To ensure that 
an kinematic constraints are resolved, only the independent set of Sij, Sj,- = 2pi ■ pj, 
can appear in this list. For example, for a four-point diagram, we can only use two 
variables of the three Mandelstam variables. 

• numeric. It is an optional variable for the basis calculation. When "numeric" is 
given and the numerical calculation in the GenerateBasis function is enabled, all the 
Grobner basis calculation is done numerically. It will speed up the computation by 
2 ~ 5 times. However, the numeric calculation has the risk of meeting kinematic 
singularities (like infrared limit and collinear limit). The numeric values should be 
rational numbers, otherwise the result depends on the floating-point tolerance inside 
the Grobner basis computation. 

• Props. This is the list for the propagator momenta. No specific order for the 
propagators is necessary. The direction of the propagator momenta is also irrelevant. 
In this version, the propagators are set to be massless. 

• RenormalizationCondition. This variable define the constraints from the renor- 
malizablity condition. Each constraint on the power of the loop momenta is expressed 
as a linear inequality. For example, when L = 3, the loop momenta are h,h,h an d 
the corresponding powers for the loop momenta are ct\ , «2 , «3 . The constraint 

ai + a 2 <6, (A.l) 

should be given as an item {{1, 1, 0}, 6} in "RenormalizationCondition" . {1,1,0} 
is the list of the coefficients of ai, 02 > «3 and 6 is the upper bound. 

The program will name the (fundamental-) scalar products (Zj • ej) as "xij". 
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A. 3 Computation and the output 

Once the input is given, all the basis determination computation is done by one command 
GenerateBasis: 

• GenerateBasis [1] . This calculates the basis analytically. 

• GenerateBasis [0] . This calculates the basis with numeric coefficients. 
The outputs are stored in the following variables, 

• ISP. This is the list for irreducible scalar products. 

• RSPSolution. This is the solutions for reducible scalar products at the unitarity 
cut. 

• CutEqnISP. This is the list for the cut equations in terms of ISPs, after all RSPs 
are eliminated. Depending on if the numeric mode is enabled, the coefficient can be 
either numeric or analytic. 

• Basis. This is the list for the terms in the basis. The output form for each term is 
(cti , . . . a m ) , where oti is the power of the i-th ISP. 

• SpuriousBasis. This is the subset of the basis which contains all spurious terms. 

• NSpuriousBasis. This is the subset of the basis which contains all non-spurious 
terms. 

• Integrand. This is the integrand after the reduction, which is an expansion over 
the integranddevel basis. cc[a±, . . . , a ni ] stands for the coefficient c Qli ,.. an for the 
term x" 1 . . . Xn" 1 , as described in (3.16). 

• Gr. This is the Grdbner basis for the ideal /, generated by cut equations in ISPs. It 
can be exported for other purposes. 

The lists "ISP" and "CutEqnISP" can be readily used for primary decomposition and 
then the dimension theory part of our algorithm, in softwares like Macaulay2 [30]. 

A. 4 Example, integrand basis for two-loop double-box diagram 

<< "/path/Basis-050712 .m" 

L=2; 

Dim=4; 

ExternalMomentaBasis= {pi , p2 , p4 } ; 

Kinematics= { pi " 2->0 , p2 " 2->0 , p4 " 2->0 , pi p2->s/2,pl p4->t/2, 
[p2 p4->- (s+t) 12, \ [Omega] l"2->-t (s+t) Is } ; 
numeric= { s->ll , t->3 } ; 

Props={ll-pl, 11, ll-pl-p2, 12-p3-p4, 12, 12 -p4, 11+12}; 
RenormalizationLoopMomenta= {{1, 0}, {0,1}, {1,1}}; 
RenormalizationPower= { 4 , 4, 6} 
GenerateBasis [1] 
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It takes about 0.95 second to generate the basis, with analytic calculation. A typical output 
is, 



Physical spacetime basis is {pi , p2 , p4 , \ [Omega] 1 } 
Number of irreducible scalar products : 4 
Irreducible Scalar Products : {xl4, x24, xl3, x21 } 

Cut equations for ISP are listed in the variable 'CutEqnISP' 
Possible renormalizable terms: 160 

The basis contains 32 terms, which are listed in the variable 'Basis' 
The explicit form of the integrand is listed in the variable ' Integrand' 
Number of spurious terms: 16 , listed in the variable ' SpuriousBasis ' 
Number of non-spurious terms: 16, listed in the variable ' NSpuriousBasis ' 
Time used: 0.955934 seconds 

The we can obtain the basis information from the variables, "CutEqnISP", "Basis", "In- 
tegrand", "SpuriousBasis" and "NSpuriousBasis". 

More examples are included in the Mathematica notebook, "examples. nb". 

B The algorithm of identifying the ISPs 

We have the following simple algorithm to find the ISPs, which is embedded in the package 
BasisDet, 

• Calculate the Grobner basis G(I') for the ideal /' generated by cut equations in terms 
of SPs, in the polynomial order of "deglex". 

• Obtain LT(G(I')), the set of the leading terms in G(I'). The linear terms in LT{G{I')) 
are the RSPs. 

It is easy to show that this algorithm gives the correct ISPs according to the definition. 

Proof. Suppose that this algorithm generates {yi, . . . ,y nR } as the list of the RSPs in the 
polynomial ordering, while {x±, . . . , x ni } as the list of the ISPs in the polynomial ordering. 
First, we can prove that y nR is a linear function of ISPs on the cut. G(I') must contain a 
linear polynomial, 



where a, hi and 7 are constants and a / 0. This polynomial cannot contain other y/s, 
because yj y y nR for j < hr. Here 'V stands for the given monomial ordering. Thus y nR 
is a linear function of the ISPs on the unitarity cut. 

Second, by induction, all yj are linear functions of ISPs on the cut. 

Third, we can prove that the ISP set is minimal. If some x% can be represented by a 
linear function of other ISPs at the cut, then 




(B.l) 




(B.2) 
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Then the leading term of this polynomial is an ISP, say x k . By the property of Grobner ba- 
sis, {LT(I')) = {LT{G{I'))). Because x k G LT(I'), x k G (LT{G{P))). Furthermore, since 
Xfc has degree one, it is generated by degree-one monomials in LT(G{I')): {y±, . . . ,y nR }- 



while 7j's are constants. This contradicts the assumption of ring structure. The ISP set is 
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